#!/usr/bin/env python 

import numpy
import random

from specks import *

s = Speck()

for i in range(2):
    s.addAtom('Al', numpy.random.normal(0, 2, 3))

beam.speck = s
beam.start()

m = Morse()

qm = QuickMin(s, m, dt = 0.05)
m.force(s)
while norm(s.f) > 10.0:
    qm.step()
    print norm(s.f)

lgv = Langevin(s, m, T = 1200.0, dt = 0.01)

totalTemp = 0.0
count = 0

while beam.running:
    lgv.step()
    count += 1
    T = 0
    for i in range(len(s)):
        T += 0.5 * s.m[i][0] * norm(s.v[i])**2    
    totalTemp += (2.0/3.0) * T / (len(s) * Kb)
    if count == 1000:
        print "Temperature averaged over 1000 timesteps:", totalTemp / count
        totalTemp = 0.0
        count = 0
    
    


